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Abstract We present an efficient, robust computational method for modeling the 
Newtonian dynamics for rotation curve analysis of thin-disk galaxies. With appropriate 
mathematical treatments, the apparent numerical difficulties associated with singulari- 
ties in computing elliptic integrals are completely removed. Using a boundary element 
discretization procedure, the governing equations are transformed into a linear algebra 
matrix equation that can be solved by straightforward Gauss elimination in one step with- 
out further iterations. The numerical code implemented according to our algorithm can 
accurately determine the surface mass density distribution in a disk galaxy from a mea- 
sured rotation curve (or vice versa). For a disk galaxy with a typical fiat rotation curve, 
our modeling results show that the surface mass density monotonically decreases from 
the galactic center toward periphery, according to Newtonian dynamics. In a large por- 
tion of the galaxy, the surface mass density follows an approximately exponential law of 
decay with respect to the galactic radial coordinate. Yet the radial scale length for the 
surface mass density seems to be generally larger than that of the measured brightness 
distribution, suggesting an increasing mass-to-light ratio with the radial distance in a disk 
galaxy. In a nondimensionalized form, our mathematical system contains a dimensionless 
parameter which we call the "galactic rotation number" that represents the gross ratio of 
centrifugal force and gravitational force. The value of this galactic rotation number is de- 
termined as part of the numerical solution. Through a systematic computational analysis, 
we have illustrated that the galactic rotation number remains within ±10% of 1.70 for a 
wide variety of rotation curves. This implies that the total mass in a disk galaxy is propor- 
tional to Vq R g , with Vq denoting the characteristic rotation velocity (such as the "flat" 
value in a typical rotation curve) and R g the radius of the galactic disk. The predicted 
total galactic mass of the Milky Way is in good agreement with the star-count data. 

Key words: galaxy: disk — galaxies: general — galaxies: kinematics and dynamics — 
galaxies: structure — methods: numerical and analytical 



1 INTRODUCTION 

Observations have shown that a galaxy is a stellar system consisting of a massive gravitationally bound 
assembly of stars, an interstellar medium of gas and cosmic dust, etc. Many (mature spiral) galaxies 
share a common structure with the visible matter distributed in a fiat thin disk, rotating about their 
center of mass in nearly circular orbits. The behavior of the stellar systems such as galaxies is believed 
to be determined by Newton's laws of motion and Newton's law of gravitation (IBinnev & Tremainel 
Il987l) . Thus, modeling the Newtonian dynamics of thin-disk galaxies is of fundamental importance 
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to our understanding of the so-called "galaxy rotation problem"-an apparent discrepancy between the 
observed rotation of galaxies an d the predictions of Newtonian dynamics as generally perceived in th e 
community of astrophysics (e.g. jFreeman & McN amara 2 0061 lRubinll2006l 120071: iBennett et al.ll2007h . 

Although scientifically well-established, the actual modeling of Newtonian dynamics, when applied 
to thin-disk galaxies, appeared in various forms in the literature with inconsistent conclusions. Without 
rigorous jus tification, some authors (e.g., lRubinll2006l 120071: IBennett et al.ll2007t ISparke & Gallagherl 
120071: iKeel 120071) tempted for simplicity to apply formulas based on Keplerian dynamics to the thin-disk 
galaxies. Theoretically, Keplerian dynamics can be derived from Newtonian dynamics as a special case 
for spherically symmetric gravitational systems such as the solar system and, therefore, is not expected 
to provide accurate descriptions for thin-disk galaxies. Hence, serious efforts were made for integrating 
the Po isson equation with mass sources distributed on a disk, as summarized by iBinnev & Tremaina 
dl 987b . The solution directly obtained from such efforts is the gravitational potential which can yield 
the gravitational force by taking its gradient. In an axisymmetric disk rotating at steady state, the gravi- 
tational force (the radial gradient of gravitational potential) is expected to equate to the centrifugal force 
due to rotation at every point. 

However, solving the disk-potential problem does not seem to be a trivial pursuit. Traditional meth- 
ods involved either treating the disk a s a flattened spheroid that consists of a serious of thin homoeoids 
each having a uniform density (e.g.. lBrandtlll960l iMestell 1 1 963 : Cuddefordl 1993b or using the sum- 
mation of modified Bessel functions f or the potential (e.g., Toomre 19631: "eemanlfl970l iNordsieckl 



Il973l ICuddefordl[l993l IConwayll200"ob . Although seemingly elegant when derived in ana lytical formu- 
las, those methods could yield closed-for m solutions only for a few special cases (e.g.. lMestellll963l 
lFreemanlfl970llBinnev & Tremainelll987h . But for determining the mass distribution in a galactic disk 
from the measured rotation curve that could have a variety of shapes, numerical integrations must be car- 
ried out and practical difficulties arise when those traditional analytical formulas are used. For example, 
the flattened spheroid approach via Abel integral and its inversion intrinsically restricts the "vertical" 
mass distribution in the disk's axial direction to that dictated by th e homoeoid structure rather than that 
from observations (e.g., according to Ivan der Kruit & Searlell 19821 the scale heights of galactic disks are 
nearly independent of radius). It is rather cumbersome to compute the surface mass density by integrat- 
ing the mass density in spheroid al shells and the "spheroid" me thods often lead to erroneous results for 
angular momentum analysis (cf . iToomrell 19631: iNordsieckl 1 973b . The Bessel function approach leads to 
an integral extending to infinity, whereas the observed rotation curve always ends at a finite distance. 
Thus, it becomes necessary to cons t ruct orbital ve l ocity beyond the observation limit b ased on various 
assumptions (e.g.. lNordsiecklfT973l lBosmalll978l: IJalocha. Bratek. & Kutscheral 12008b . Moreover, the 
derivative of rotation velocity usually appearing in the Bessel function formulation for computing mass 
density tends to introduce significant errors in practical applications. 

In general, the fundamental solution to th e Poisson equation (tha t governs the gravitation poten- 
tial) is called Green's function dArfkeniri985t ICohl & Tohlind[T999b . The potential from arbitrarily 
distributed sources can be obtained by integrating the Green's function-serving as the integral kernel- 
multiplied by the source density throughout the region where the sources are located. Thus, considering 
the gravitational potential in terms of Green's function is the most direct approach for realistic modeling 
the g alactic rotation dynamics (e.g., lEckhardt & Pestaflall2002l : iPierens & Hurell2004t Ih ure & Pierensl 
12005b . For sources distributed axisymmetrically on a thin dis k, the Green's function can be expressed 
in terms of the complete elliptic integral of the first kind (e.g. jBinnev & Tremaindfl987b . Because the 
dynamics of thin-disk galactic rotation is typically described along the midplane (z = 0) with the mass 
distribution being symmetric about the disk midplane and about its central axis, the radial gradient of 
potential in the midplane must be evaluated. The elliptic integrals of the first kind and second kind that 
appear in the radial gradient of potential can become mathematically singular at the midplane (when 
z = 0) where the radius of the source approaches that of the observation p oint. Such singularities have 
been considered "inco nvenient from the po i nt of v iew of numerical work" bv lBinnev & Tremainel(ll987l) 
and "bothersome" by lEckhardt & Pestanal (|2()()2). Methods were suggested to circumvent such singu- 
lar itie^>_bj^valuatingjh^ vertical distance z slightly away from z = 
(cf lBinnev & Tremainel 19871 lEckhardt & Pestafiall2002h . which seem to be somewhat ad hoc by nature 
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and lack of desirable mathematical elegance. On the other hand, it is the axisymmetric mass distribution 
within an idealized rota ting infinitesm ally thin disk that has often been of practical interest especially for 
rotation curve analysis dToomrel 19631) . Therefore, the efforts of effectively dealing with the singularities 
arising from elliptic integrals has been continuously made f or robust and accurate computations of the 
disk galaxy rotation problem (eve n up to recent years, e.g., lEckhardt & Pestaflall2002t iPierens & Hurel 
l2004ilHure & Pierenj2005ll2009l) . 

In the present work, we derive a numerical model for computing the Newtonian dynamics of thin- 
disk galactic rotation that allows the mass to be distributed even in an infinitesmally thin region around 
the midplane of the disk with the governing equation being considered strictly along the midplane 
(z = 0) and the singularities from elliptic integrals treated rigorously based on the concept of the 
mathematical limit. To enable dealing with arbitrary forms of rotation curves and mass density distri- 
butions, we adopt the techniques developed w ith boundary element method (cf. lSladek & Sladekll 19981 : 
Idrav 111 9981: [Sutradhar. Paulino & Gravll2008l) for solving integral equations using compactly supported 
basis functions instead of that extending to infinity like Bessel functions, as detailed in Section 2. Hence 
the finite physical problem domain for disks of finite sizes can be conveniently considered, without the 
need of speculated rotation curve beyond the "cut-off" radius and evaluation of the derivative of rota- 
tion velocity. By nondimensionalizing the governing equations, a dimensionless parameter which we 
call the "galactic rotation number" appears in the force balance (or centrifugal-equilibrium) equation, 
representing the gross ratio of centrifugal force and gravitational force. We show that together with a 
constraint equation for mass conservation, the value of this galactic rotation number can be determined 
as part of the numerical solution, with computational examples presented in Section 3. With a known 
value of the galactic rotation number, the total galactic mass can be determined from measured galactic 
radius and characteristic rotation velocity, as shown in Section 4 wherefrom important physical insights 
are discussed. 



2 MATHEMATICAL FORMULATION AND COMPUTATIONAL TECHNIQUES 

For convenience of mathematical treatment, we represent a rotating galaxy by a self-gravitating contin- 
uum of axisymmetrically distributed mass in a circular disk with an edge at finite radius R g , as shown in 
Fig-El This kind of continuum representation is typically valid when the distributed masses are viewed 
on a scale that is small compare to the size of the galaxy, but large compared to the mean distance be- 
tween stars. Without loss of generality, we consider the thin disk having a uniform thickness (h) with 
a variable mass density (p) as a function of radial coordinate (r). Because we consider the situation of 
thin disk, the vertical distribution of mass (in the z-direction) is expected to contribute inconsequen- 
tial dynamical effect especially as the disk thickness becomes infinitesmal. In mathematical terms, the 
meaningful variable here is actually the surface mass density er(r) = p(r) h. Whether to consider the 
surface mass density er(r) or the bulk mass density p(r) in the mathematical equations is really a matter 
of taste, since they can easily be converted to each other using a constant factor h by our definition. In 
the present work, we use the bulk density p(r) for its consistency with the direct physical perception of 
a thin disk with a nonzero thickness h. 

Physically, the stars in a galaxy must rotate about the galactic center to maintain the disk-shape 
mass distribution. Without the centrifugal effect due to rotation, the stars wo uld collapse into the galacti c 
center as a result of the gravitational field among themselves. According to lBinnev & Tremainel(ll987l) . 
it is also reasonable to assume the galaxy is in an approximately steady state with the gravitational force 
and centrifugal force balancing each other, in view of the fact that most disk stars have completed a 
large number of revolutions. 



2.1 Governing Equations 

Instead of following the traditional approach by first solving gravitational potential from the Poisson 
equation, we derive the governing equation directly from the consideration of force balance. Here, the 
force density on a test mass at the point of observation (r, = 0) generated by the gravitational attraction 
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Fig. 1 Definition sketch of the thin-disk model considered in the present work. The mass 
is assumed to distribute axisymmetrically in the circular disk of uniform thickness h with a 
variable density as a function of radial coordinate r (but independent of the polar angle 9). 



due to the summation (or integration) of a distributed mass density p(f) at position described by the 
variables of integration (f , 9) is expressed as an integral over the entire disk, with the distance between 
(r, 9 = 0) and (f, 6) given by (f 2 + r 2 — 2f r cos 9) 1 / 2 and the vector projection given by (r cos 6 — r). 
Thus, the equation for gravitational force to balance the centrifugal force at each and every point in a 
thin disk can be written as (according to Newton's laws) 
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{fcosd — r)d6 
o (f 2 + r 2 — 2frcos6) 3 / 2 



p{f)hfdf 



V(r) 



(1) 



where all the variables are made dimensionless by measuring lengths (e.g., r, f, h) in units of the 
outermost galactic radius R g , disk mass density p in units of M g /R g with M g denoting the total galactic 
mass, and rotation velocities V(r) in units of the a characteristic galactic rotational velocity Vq (usually 
defined according to problem of interest, e.g., such as the maximum velocity corresponding to the fiat 
part of a rotation curve). The disk thickness h is assumed to be constant and small in comparison with the 
galactic radius R g . Our results for surface mass density p(r) h are expected to be insensitive to the exact 
value of this ratio as long as it is small. There is no difference in terms of physical meaning between 
the notations (r, 9) and (f, 9); but mathematically the former denotes the independent variables in the 
integral equation (Q~|i whereas the latter the variables of integration. The gravitational force represented 
as the summation of a series of concentric rings is described by the first (double integral) term while the 
centrifugal force by the second term in ([T). 

Our process of nondimensionalization of the force-balance equation yields a dimensionless param- 
eter, which we call the "galactic rotation number" A, as given by 



Vq 2 R 9 

Mg G 



(2) 



Modeling Newtonian Dynamics of Thin-Disk Galaxies 



5 



where G (= 6.67 x 10~ n (m 3 kg -1 s -2 )) denotes the gravitational constant, R g is the outermost galactic 
radius, and Vq is the characteristic velocity (which is equated here to the maximum asymptotic rotational 
velocity). This galactic rotation number A simply indicates the relative importance of centrifugal force 
versus gravitational force. 

Equation ([TJ can either be used to determine the surface mass density p(r) h from a given rotation 
curve V(r) or vice versa. But when both p(r) and A are unknown, another independent equation is 
needed to have a well-posed mathematical problem. According to the law of conservation of mass, the 
total mass of the galaxy M g should be constant satisfying the constraint 

2tt / p(f)hfdf = 1. (3) 
Jo 

This constraint can be used for determining the value of galactic rotation number A while (Q]) for p(r). 
Equations ([Hi-tO can in principle be used to determine the mass density distribution p(r) in the disk, 
the galactic rotation number A, and the total galactic mass M g , all from measured values of V(r), R g , 
Vo, and h. On the other hand, if p(r) and h (or p(r) h) are known, V(r) can of course be determined 
from ((TJ. 

Moreover, it is known that the integral with respect to 9 in ([TJ can be written as 

r2w 



(f cos6* — r)d6 



K(m) 



E(m) 



r(f — r) r(f + r) 



lo (f 2 + r 2 — 2fr cos 6*) 3 / 2 
where K(m) and E(m) denote the complete elliptic integrals of the first kind and second kind, with 

4fr 



(4) 



(f + r) 2 



Thus, (JJi can be expressed in a simpler form 

° l r E{m) K{m) 



p(r)hrdr + -AV{r) 2 



0. 



(5) 



(6) 



which is more suitable for the boundary element type of numerical implementation (with the double 
integral converted to a single integral). 



2.2 Computational techniques 



Following a standard boundar y element approach (e.g., ISladek &Sladekl Il998t iGravl 119981: 
ISutradhar. Paulino & Gravll2008l) . the governing equations d6) and (O can be discretized by dividing 
the one-dimensional problem domain < r < 1 into a finite number of line segments called (linear) 
elements. Each element covers a subdomain confined by two end nodes, e.g., element n corresponds to 
the subdomain [r n , r n+ i], where r n and r n+ i are nodal values of r at nodes n and n+1, respectively. On 
each element, which is mapped onto a unit line segment [0, 1] in the ^-domain (i.e., the computational 
domain), p is expressed in terms of the linear basis functions as 

P(0 = Pn(l - + Pn+lL 0<e<l, (7) 

where p n and p n +i are nodal values of p at nodes n and n + 1, respectively. Similarly, the radial coordi- 
nate f on each element is also expressed in terms of the linear basis functions by so-called isoparameteric 
mapping: 

r(0=r n (l-0+r n +iZ, 0<£<1. (8) 
If the rotation curve V(r) is given (from observational measurements), the TV nodal values of p n = 
p{r n ) are determined by solving N independent residual equations over N — 1 element obtained from 
the collocation procedure, i.e., 



JV-l 

E 

n=l 1 



E(mi) K(rrii) 



df 1 

Pmm^dt + -AV{ ri ) 2 = , i = 1,2, N , 



(9) 
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with 

™«®= r^ ( i ) 'V (10) 

where p(£) = p„ (1 — £) + p n +x £■ The value of A can be solved by the addition of the constraint equation 

jv-i r i 

civ 

p{£)M{£)-dti-l = Q. (11) 



n=l 

Thus, we have N + 1 independent equations for determining N + 1 unknowns. The mathematical 
problem is well-posed. The set of linear equations comprising (0 and dTTb for N + 1 unknowns (i.e., N 
nodal values of p n and A), once computed with appropriate treatments of the mathematical singularities 
shown in Appendix A, can be transformed into a matrix form using the Newton-Raphson method and 
then solved with a standard matrix solver, e.g., by Gauss elimination in one step without further iterations 
dPress et al.ll 19881) . 

3 COMPUTATIONAL EXAMPLES 

As we mentioned before, equations 01 and ( fTTT i can be used to either solve for p(r) and A from a given 
rotation curve V(r) or determine the rotation curve V(r) from a given surface mass density distribution 
a{r) = p(r) h. Usually, solving for p(r) from a given rotation curve V(r) requires computation of a 
linear algebra matrix problem whereas determining V(r) from a given p(r) only involves a straightfor- 
ward integration. But in a spiral galaxy it is the rotation curve that can be measured with considerable 
accuracy; therefore, the observed rotation curve has been r egarded to provide the most reliable m eans for 
determining the distribution of gravitating matter therein (lToomrelll963l:ISofue & Rubinll200ll) . Hence, 
we first consider examples of solving for p(r) and A from a given V(r). 



3.1 Mass distribution for rotation curve of typical shape 

To obtain numerical solutions, the value of (constant) disk thickness h must be provided; we assume h = 
0.01, which is typical of disk galaxies like the Milky Way. For computational efficiency, we distribute 
more nodes in the regions (e.g., near the galactic center and disk edge) where p has a greater gradient 
of variations. The typical number of nonuniformly distributed nodes N used in the computation is 1001 
with which we found for most cases to be sufficient for obtaining a smooth curve of p versus r and 
discretization-insensitive values of galactic rotation number A. When numerically integrating element- 
by-element in (0 and dTTb , we use ordinary 6-point Gausian quadrature for integrals with respect to 
< £ < 1. The two-dimensional integrals JA.4-b on a singular element are calculated numerically by 
ordinary 6 x 6-point Gausian quadrature on a unit square with < rj < 1 and < £ < 1. 

The measurements of galactic rotation curve of mature spiral galaxies reveal that the rotation 
velocity V(r) typically rises linearly from the galactic center in a small core and then bend down 
to reach an approximately constant value extending to the galactic perip hery dRubin & Fordl [l970t 
iRoberts & Whitehurstl [19751: |Bosmalll978l; iRubin, Ford. & Thonnardlll98ol) . These basic features may 
be mathematically idealized as 

V(r) = 1 - e- r/R ' , (12) 

where the dimensionless orbital velocity V(r) is measured in units of the characteristic velocity Vb 
defined as the maximum orbital velocity, and the parameter R c can serve as the scale of the "core" of 
a galaxy. As shown in Fig [2] close to the galactic center when r/R c is small, we have V(r) ~ r/R c 
describing a linearly rising rotation velocity (by virtue of the Taylor expansion of e~ r / Rc ). The initial 
slope of this rising rotation velocity is given by l/i? c . Thus, larger value of R c leads to a more gradual 
rise of the rotation velocity and a shrinking "fiat" part of rotation curve which seems to disappear for 
R c > 0.2. 

Corresponding to the rotation curves in Fig |2] as described by (fT2l . the computed mass density 
distributions in galactic disk are shown in Fig. [3] For R c < 0.02, the curves of p versus r approach an 
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Fig. 2 Nondimensionalized orbital velocity profiles V(r) according to mathematically ide- 
alized description (US for R c = 0.01, 0.02, 0.05, 0.1, and 0.2. 



asymptotic one for the most part except in a tiny region around galactic center where the peak density 
value at r = still increases with further decreasing R c . In other words, the mass density tends to 
decrease rapidly from the galactic center (with a slope becoming steeper for a tigher galactic core with 
a smaller R c ). However, beyond r = R c , the mass density decrease more gradually towards the galactic 
periphery until reaching the galactic edge where it takes a sharp drop. Outside the galactic core (r > R c ), 
only for R c > 0.1 do changes in mass density distribution and the value of A become noticeable with 
varying R c . Noteworthy here is that the computed values of galactic rotation number A for R c < 0.15 
are within a small interval [1.5708, 1.6422] despite orders of magnitude of R c variation. It appears that 
as R c — > the value of A approaches a limit at ~ 1.5708. For example, the computed results show 
that A = 1.57085 and 1.57080 for R c = 0.005 and 0.001, respectively. But the increase of A with R c 
becomes more significant for R c > 0.15, as illustrated by the computed results at R c = 0.2 and 0.3 
yielding A = 1.7098 and 1.9224, respectively. 

At the limit of R c — > 0, the (idealized) rotation curve as described by (Q~2) approaches a completely 
flat one V(r) = 1 (except in the infinitesmal n eighborhood o f r = 0). The solution at this limit should 
approach that of the well-known Mestel's disk dMeste ill 19631) given by 



p{r) 



A 
2-nhr 



1 - -sin -1 (r) 

7T 



(13) 
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Fig. 3 The distributions of mass density p(r) computed for R c = 0.01, 0.02, 0.05, 0.1, and 
0.2 with A = 1.5710, 1.5719, 1.5777, 1.5999, and 1.7098 determined as part of the numerical 
solutions. 



in a dimensionless form consistent with the nomenclature of the present work. Here, according to (|3) 
the galactic rotation number A can be determined by 

A=— j = - = 1.5707963. (14) 

^[l-fsin-^r)] df 2 

As a test, we can substitute p(r) given by $13[ into (O and ( fTTT i and compute with our code for 
numerical integrations to determine V(r) and A. With the first node at r = being ignored to avoid the 
numerical difficulties with the singularity of p in ( TT3l l. we can indeed obtain a flat V(r) = 1 throughout 
the entire interval (0, 1] (except in an infinitesmal neighborhood around r = 0) and A = 1.57081. The 
computed curve of p versus r corresponding to (fT3) with A = 1.57081 overlaps that of R c = 0.01 in 
Fig [3] (except in the infinitesmal neighborhood of r = 0), as expected. This exercise demonstrates our 
code capability for determing the rotation curve from a given disk mass distribution, and also in a way 
verifies the correctness of our computational code implementation. Since most Sb galaxies-intermediate 
type of spiral galaxies-have rotation curves typically with a very steep rise in a small central core region, 
the mass density distribution in those Sb galaxies (including the Milky Way) is expected to be reasonably 
approximated by that of the Mestel disk dT3"V But for less massive Sc galaxies having more gradual rise 
rotation curves, their mass density distributions can deviate noticeably from that of the Mestel disk 
especially toward the galactic center, as shown in Fig|3]for those with R c > 0.02. 
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3.2 Rotation curve for given mass distribution 

As demonstrated in Section 3.1, numerically computing the integration in (0 for a given p[f) as that of 
Mestel's disk can produce a completely flat rotation curve. Actually, rotation curves similar to those in 
Fig [2] can also be produced by a combination of the Freeman exponential disk and Mestel disk. Here, 
the Freeman disk has a surface mass density proportional to e~ r l Rd with R c i denoting a scale length for 
the exponental disk dFreemanlll970l) . But the Freeman exponential disk alone is known not to be able 
to produce a rotation curve with considera ble fiat portion as often being observed in disk galaxies (e.g., 
lFreemanlll970H B innev & Tremainelll987l) . The case of V(r) for the Freeman exponential disk can also 
be computed with our code, as a check; the result showed excellent agreement with that of Freeman's 
analytical formula. If we use the Freeman disk for describing the galactic core having a rising rotation 
velocity and Mestel disk for the outer flat part, there is a good chance to obtain rotation curves of 
typically observed shapes. For example, we can simply construct a mass density model (which we call 
the Freeman-Mestel model) as 



p(r) = 



Po e 



-r/R d 



< r < R c 



^[l-fsin-V)] , R c <r<l 



(15) 



where 



and 




Rd = < 



Po 



Wl - i?2[l - 2sin- 1 (i? c )/ 7 r] 



l--sm- 1 {R c ) 



2-KhR r e- R ol R <i 



so that both p and dp/dr are continuous at the connecting point r = R c . Moreover, the mass conserva- 
tion constraint (fTTT i can be used to determine the value of galactic rotation number as 



.4 



JV-l 



df 



2n E J o p*(Ohm^ 



(16) 



where p* comes from that given by (T5[ by setting A = 1 . 

Although R c here also serves as a scaling parameter for the galactic core, having a similar physical 
meaning as R c in (fl~2l . the value of R c does not have any mathematical relationship with that of R c . 
For example, at R c = 0.05 ( f]~5b and ( fl~6] l yield V(r) and p(r) in Figs[4]and[5]noticeably different from 
those in Figs [2] and [3] For smaller values of R c , the differences between p(r) given by the Freeman- 
Mestel model and that in Fig[3]at the same values of R c are less visually discernable. But the value of A 
determined by the Freeman-Mestel model can still be slightly different. For example, at R c = R c = 0.01 
( [Tol l yields A = 1.5777 whereas that computed in Section 3.1 is A = 1.5710. It seems for a given value 
of R c = R c the rotation curve of the Freeman-Mestel model has a greater slope for the rising velocity 
in galactic core but a somewhat less flat velocity outside the core, as shown in Fig [4] Such a numerical 
difference tends to diminish with diminshing R c , e.g., we have A = 1.57147, 1.57084, and 1.57081 for 
R c = 10~ 3 , 10~ 4 , and 10~ 5 , respectively. As expected, A —> 1.57080 as that for the Mestel disk given 
in 03} at the limit of R c -> 0. 

What we try to illustrate here is that for obtaining rotation curves with basic observed feat ures, 
a simple analyt ical mass density model as constructed by combination of those of iMesteil d 19631) and 
iFreemanI d 19701) in (TT3T > seems to be quite reasonable and convenient. In terms of computational efforts, 
it is usually much easier and faster to compute the rotation velocity V(r) from a given mass density 
distribution p(r) than vice versa. This is because that computing V(r) for a known p(r) does not need 
to solve the matrix problem. However, there has not been reliable means for directly measuring the 
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Fig. 4 The rotation velocity V(r) determined with p(r) given by (25) for the Freeman-Mestel 
model at R c = 0.05, compared with that in Fig[2]for R c = 0.05. 



mass distribution in a disk galaxy. The mass distribution derived from measured luminosity must rely on 
assumed mass-to-luminosity ratio, with the validity of which being a subject of debate. Thus, accurately 
measured rotation curves remain as the most reliable basis for determining the distribution of mass 
in dis k galaxies, providing f undamental information for understanding the stellar dynamics in galactic 
disks dSofue & Rubinll200ll) . 

3.3 Analysis of measured rotation curves of arbitrary shapes 

For rotation curves with "idealized" shapes expressed in terms of simple mathematical functions like 
that in dT2b . we have shown that the numerically computed mass density distribution p(r) approaches 
that of Mestel's disk dot when the galactic core is small, e.g., for R c < 0.02. But some measured 
rotation curves can vary significantly from those described by simple mathematical functions or those 
produced by conveniently constructed mass density functions like with the Freeman-Mestel model (TT3T >. 

To determine the mass density distribution according to Newtonian dynamics from a measured ro- 
tation curve of arbitrary shape, our computational scheme based on sound mathematical foundation as 
presented in Section 2 (as well as Appendix A) can become a generally applicable and flexible tool for 
many practical applications. As an example, here in Figs [6] and [7] we show our computed mass density 
distributions for a few actually measured galactic rotation curves with noticeably different characteris- 
tics. 
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Fig. 5 The distribution of mass density p(r) given by (25) for the Freeman-Mestel model at 
R c = 0.05 with A = 1.6060, compared with that in Fig[3]for R c = 0.05 with A = 1.5777. 



The measured rotation curve for the Milky Way in Fig|6]seems to be just a few bumps and wiggles 
superposed on that in Fig[2]for R c = 0.01. Therefore, it is no surprise to see that the corresponding mass 
density curve for the Milky Way in Fig Q also exhibits a few bumps and wiggles around that in Fig [3] 
for R c = 0.01. Similarly, the measured rotation curve for NGC 3198 in Fig appears to be just that in 
Fig|2]for R c = 0.05 with some small perturbations, and so does the computed NGC 3198 mass density 
in Fig|7]compared with that for R c = 0.05 in Fig[3] But the rotation curve for NGC 2708 in Fig |6]differs 
significantly from those of typical shapes in Fig [2] The computed mass density distribution for NGC 
2708 in Fig [7] shows noticeably different features from those in Fig [3] The sharp rise of mass density 
forward galactic center corresponds to a fast dropping of rotation velocity, as required for the force 
balance in Newtonian dynamics. The gradual increase in the rotation velocity in the middle section 
(0.1,0.7) of NGC 2708 leads to a slow decreasing local mass density. Then a slight reduction of the 
rotation velocity toward the galactic periphery is responsible for a faster decrease of local mass density 
in the outer region r > 0.7 than those for flat rotation curves in Fig |7]for NGC 2708. 

Despite the differences in rotation curves in Fig [6] the computed values of galactic rotation number 
A for these three galaxies are quite close within a few percents, namely, A — 1.564, 1.619, and 1.644, 
respectively for the Milky Way, NGC 3198, and NGC 2708. This is consistent with that shown in Fig[3] 
for a wide range of R c . Thus, we may reasonably conclude that for most disk galaxies, the value of 
galatic rotation number is expected to be within ±10% of A = 1.70, with smaller A for the galaxies 
having high-density core and small R c and larger A for those having more gradual rise in the rotation 
curve with larger R c . 
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Fig. 6 The rotation curves V(r) of Milky Way with V = 2.2 x 10 5 (m s _1 ) and R g = 
4.73 x 10 20 (m), NGC 3198 with V = 1.5 x 10 5 (m s" 1 ) and R g = 9.24 x 10 20 (m), and 
NGC 2708 with V = 2.3 x 10 5 (m s" 1 ) and R g = 1.42 x 10 20 (m). 



Although we only computed examples with a few representative rotation curves, such as those 
described by the idealized formula (fTZt with several values of R c and those actually measured with 
different characteristics, we believe the cases examined here actually cover a wide enough range of 
observational measurements that our results can offer general physical insights. Cases with rotation 
curves falling either close to or in between those illustrated here are not expected to differ considerably 
from our present findings. 

4 DISCUSSION 

The problem of determining the mass distribution in a thin axisymmetric disk from observed circular 
velocities has been investigated by many authors over the past fifty years, through various mathemat- 
ical approaches. Yet satisfactory method for accurate computation is still lacking, despite the galactic 
rotation model has been simplified as much as possible for concisely describing only the most essential 
features. The main obstacle here appears to have been due to the mathematical singularities in the el- 
liptic integrals that are apparently difficult to handle. Here in this work, we present an efficient, robust 
computational method with appropriate mathematical treatments such that the apparent difficulties as- 
sociated with the singularities are completely removed. Thus, we are enabled to systematically analyze 
the basic features in a rotating disk galaxy, with properly nondimentionalized mathematical formula- 
tions. Further refinement of the present galactic rotation model may provide description of some o f the 
fine details such as the spiral arm structure with non-axisymmetic motion dKoda & Wa da 2002]), gas 
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Fig. 7 The computed mass density distributions p(r) from given rotation curves in Fig |6]for 
Milky Way, NGC 3198, and NGC 2708, with the values of A determined as 1.564, 1.619, and 
1.644, respectively. 



pressure effect in the central core dDalcanton & Stilpll2010l ). disk thickness effect dCasertanoll 1983b . etc. 
But those fine details should not alter the basic features significantly, at least in a gross qualitative sense. 
Our results in Fig 7 show that the general shape of the mass density distribution remains quite similar 
for rotation curves of drastically different appearances. The value of the galactic rotation number A does 
not change more than ±10% for a variety of rotation curves, indicating that the gross balance between 
the centrifugal force and gravitational force in a disk galaxy is usually insensitive to fine details. 



4.1 Total mass in galactic disk 

In the dimensionless form as presented here, our mathematical system contains a dimensionless param- 
eter called galactic rotation number A. This galactice rotation number, with its value determined as part 
of the computational solution, can provide a unique insight into the dynamical system of rotating galaxy. 
From the knowledge of Vq and R g from measured rotation curves, we can determine the value of total 
mass M g based on computed value of A (cf. (O) as 



AG 



M a = ~a~7^~ ■ (17) 



According to the rotation curve of the Milky Way in Fig 6, we have the galactic rotation number 
A = 1.564. Then, from the measured values Vq = 2.2 x 10 5 (m s _1 ) and R g = 5 x 10 4 (light-years) 
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= 4.73 x 10 20 (m) (which is about 15.3 kpc where 1 kpc = 3.086 x 10 19 m), O yields 

M g = 2.19 x 10 41 (kg) = 1.10 x 10 11 (solar-mass) . (18) 

(Here, 1 solar-mass = 1. 98892 x 10 30 kg.) This val ue is in very good agreement with the Milky Way 
star counts of 100 billion dSparke & Gallagherll2007l) . 

Another example in Figs 6 an d 7 is the galaxy NGC 3198, with Vq = 1.5 x 10 5 (m s _1 ) and 
R g = 30 (kpc) = 9.24 x 10 20 (m) (lBegemanlll987. 1989h . Using the computed A = 1.619, we obtain 
M g = 1.925 x 10 41 (kg) = 9.68 x 10 10 (solar-mass). 

For a small disk galaxy NGC 6822, we have a rotation curve similar to th at described by R r ~ 0.3 
in O, with V = 6.0 x 10 4 (m s" 1 ) and R g = 5 (kpc) = 1.54 x 10 20 (m) dWeldrake et al.ll2003l) . If 
we take A = 1.92 for R c = 0.3, (El) yields M g = 4.33 x 10 39 (kg) = 2.18 x 10 9 (solar-mass). 

Because the value of A does not vary much for a large range of rotation curves with va rious shapes 
(see, e.g., Figs [6] and [7]), what (JT7J implies is that M g oc V 2 R g as what iBosmal d 1 978b found from 
evaluating mass versus size in a large number of observed disk galaxies. For a fixed value of Vq, M g oc 
R g . Therefore, a disk galaxy cannot physically extend indefinitely in size, for M g to remain finite. In 
other words, there must be an edge of the galactic disk at a finite radius R g , where the mass density 
precipitously diminishes. Normally, one would define R g as the radial distance where the "luminous", 
"visible", or "detectable" signal for rotating matter ends. With the advance in measurement technology 
using different emission lines, the dete ctable rotating matter ( in the form of gas) seems to extend further 
out from the optically visible disk (cf. ISofue & Rubinll200ll) . Thus, the value of R g may change with 
the evolving astronomical observation technology. Wherever the true R g is located, it must correspond 
to an abruptly steep decrease of mass density whereas the mass density variation within R g is expected 
to be smooth, according to our Newtonian dynamics model for thin-disk galaxies with typical rotation 
curves. It should be noted that although for a given rotation curve with fixed Vq the total mass M g of the 
galactic disk increases linearly with R g , the dimensional value of surface mass density should generally 
decrease with R g according to 1/R g because it scales as M g /R g . 

As an interesting exercise, we may take (fl~3~b for the convenience in estimating the surface mass 
density a(r) = p(r) h around the Sun in the Milky Way when R g increases. Then, we obtain a(r sun ) = 
0.3106, 0.7954, and 1.7532 for r sun = 0.5229, 0.2614, and 0.1307, respectively for R g = 15.3, 30.6, 
and 61.2 (kpc) assuming the Sun is located at r sun R g = 8 (kpc) from the galactic center. Based on 
the value given by ( fT8b , we have the dimensional surface mass density a(r sun ) M g /R 2 w 146 (solar- 
mass pc -2 ) for R g = 15.3 (kpc). If R g for the flat rotation curve were found to be at 30.6 or 61.2 
(kpc), the dimensional surface mass density would become 187 or 206 (solar-mass pc~ 2 ), varying much 
less dramatically than the value of R g . This phenomenon is a consequence of the 1/r part of (fT3l . 
which becomes more dominant for smaller values of r. In fact, if the surface mass density er(r) were 
strictly to follow a distribution oc 1/r, the dimensional surface mass density for a given dimensional 
radial coordinate rR g would remain constant because the value of A changes little if at all. Thus, as 
Rg extends further out, the value of dimensional surface mass density in the neighborhood of Sun is 
expected to become almost independent of the value of R g . 

4.2 Computed mass density versus observed surface brightness 

Observations of disk galaxies suggest that the surface brightness-the total stellar l uminosity emit- 
ted per unit area of the disk-is approximately an exponential function of radius dFreemanl Il970l 
iBinnev & Tremainelll987l) . This exponential approximation seems to b e especially go od for the outer 
part of disk galaxies where the inner bulge component diminishes (e.g.. lFreemanlll970h . Our computed 
mass density distributions in Fig [3] according to typical flat rotation curves indeed show nearly straight- 
line shape in the log-linear plots corresponding to approximately exponential function for a large portion 
of the problem domain, e.g., in the interval (0.2, 0.9). In fact, the least-square fit of our computed In p 
versus r for the case of R c = 0.01 (cf. FigO to a linear function for 0.2 < r < 0.9 yields 



In p= 5.2614 - 3.4377 r, 



(19) 
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with the correlation coefficient "i? 2 " being 0.9968 suggesting that the portion of mass density in 
(0.2, 0.9) can indeed be well described by an exponential function p = po e~ r l Rd with po = 192.75 
and Rd = 0.2909. If the same least-square fitting were done for 0.1 < r < 0.9, we would have 
po = 238.41 and R d = 0.2668 but with a slightly reduced correlation coefficient R 2 = 0.9870, which 
still indicates a good approximation with the exponential function. However, the dimensional "radial 
scale length" R d R g for the Milky Way would be ~ 4.5 (or 4.1) (kpc) according to R d = 0.2909 (or 
0.2668) assuming R g = 15.3 (kpc). This is larger than the rad ial scale length 2.5 (kpc) from fitting the 
brightness measurement data reported by iFreudenreichl (ll 99 8l) . For NGC 3198 with R g = 30 (kpc), we 
would have R d R g = 8. 73 (or 8.00) (kpc), ag ain larger than the radial scale length of 2.63 (kpc) for 
the luminosity profile (cf. lBegemanll987. 19891) . So, our computed results suggest that the surface mass 
density decreases toward the galactic periphery at a slower rate than that of the luminosity density. In 
other words, the mass-to-light ratio in a disk galaxy is not a constant; it generally increases with the 
radial distance from the galactic center as indicat ed by our anal ysis for the exponential portion of mass 
density distribution (which was also suggested bv lBosmd[l978h . 

But it is known that the constructed mass density distribution in terms of a single exponential func- 
tion cannot generate an observed flat rotation curve (lFreemanlll970t iBinnev & Tremainelll987l) . The 
sharp increase of the mass density near the galactic center that drastically deviates the exponential de- 
scription for 0.1 < r < 0.9 or 0.2 < r < 0.9 seems to play an important role for keeping the rotation 
curve flat forward the galactic center up to the edge of the core. In reality, most disk galaxies also have a 
central bulge with apparently high concentration of stars. Our pure disk model does not explicitly treat 
the bulge as a separate object; instead, the gravitational effect of the bulge is lumped in the rotating disk. 
Thus, our computed mass density should be regarded as a combination of that from the pure disk and 
the effective bulge represented in the disk form. This sharp increase of the disk mass density near the 
galactic center can be considered as an account for the highly concentrated mass in the central bulge. 
Actually, it may not be impossible to extend the formulation in Section 3.2 for a mass density distri- 
bution to include a summation (or expansion) of several exponential terms with different radial scales 
lengths, for matching an observed rotation curve with more complicated shape. Yet, the most straight- 
forward method for determining the mass density distribution for a given rotation curve (of arbitrary 
shape) is by numerically solving the linear algebra matrix equation derived based on sound mathemati- 
cal ground for disk galaxies of finite size as presented in Section 2 and demonstrated in Section 3.1 and 
Section 3.3. 

4.3 Inaccuracy of Keplerian dynamics for disk galaxies 

Enchanted by its simplicity, the Keplerian dynamics was applied by several authors in description of 
the disk galaxy behavior without seriously inquiring its validity and accuracy. To clarify some of the 
problems in such an over-simplification, here we present a quantitative analysis of the fundamental 
differences between the Keplerian dynamics and Newtonian dynamics especially when applied to disk 
galaxies. 

From analyzing the orbits of planets around the Sun, Kepler empirically discovered laws for planet 
motion in the solar system. It was Newton who mathematically showed that Kepler's laws are actually 
consequences of Newton's laws of motion and universcal law of gravitation. The so-called Keplerian 
dynamics can be derived from Newton's theorems for the gravitational potential of any spherically 
symmetric mass distribution. In considering the balance between gravitational force from the distributed 
mass in a galaxy and centrifugal force due to rotation, applying Keplerian dynamics would lead to an 
equation as 



which is apparently quite different from (fl3 as rigorously derived for the thin-disk galaxies. However, the 
force balance equation based on Keplerian dynamics ( f20b looks much simpler than that of Newtonian 
dynamics (Q]). If justifiable in a quantitative sense, it may be conveniently used as a reasonable approxi- 




(20) 
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mation to the more involved rigorous computations. To provide a quantitative comparison, we herewith 
examine a few basic mathematical features of ( l20b to illustrate whether the Keplerian dynamics can be 
practically used as a reasonable approximation to the Newtonian dynamics (Q]i for disk galaxies. 

For a given rotation curve with the orbital velocity V(r) described by ( fT2l . an analytical solution to 
( f20l > for p(r) can be obtained as 



p(r) = 



A 
2^~h 



- ( 1 - 2e- r l R - 
r 



,-2r/R c 



2 

R~c 



e -r/R c _ e -2r/R c 



(21) 



Thus, ( f2TT > describes a mass density approaching 3Ar/(2ir h i? 2 ) — > as r — > with a positive slope 
for small r yet approaching Aj (2tt h r) as r — > 1 (because e~ 1 ' R " can be negligibly small for small R c , 
e.g., e - x ' Rc = 4.54~ 5 , 2.06 x 10" 9 , and 1.93 x 10~ 22 for R c = 0.1, 0.05, and 0.02, respectively). The 
mass density distribution of (|2TT > does not monotonically decrease with r as that shown in Fig[3j instead, 
it is zero at the galactic center and increases for small r according to a slope oc 1/i? 2 (which can be 
large for small R c ) until reaching a peak value, then decreases in a form oc 1/r towards the galactic 
periphery r = 1 without the precipitous drop seen in Fig [3] 
Substituting (EB to (O yields 



.4 



1 



1 - 2e~ r / R - 



-2r/R c 



(22) 



which leads to A as 1 for small R c , quite different from 1.57 when R c — > as obtained in Section 3.1. 
Hence using the Keplerian dynamics to describe the disk galaxies can be misleading, because not only 
the results differ quantitatively but also qualitatively from that based on rigorous computations. 

On the other hand, if we assume the mass density distribution is known, e.g., as that given by ( [T3l >. 
(f20l > leads to 



V{rf 



1 



1 sin (f) 

7T 



df = l- 



sin (r) — 



l - vi^T 5 



Instead of a completely flat rotation curve, the Mestel's disk mass density distribution with Keplerian 
dynamics would yield orbital velocity V(r) that monotonically decreases with r, having V(0) = 1 
and V(l) = 0.7979. Therefore, a mass density distribution corresponding to a flat rotation curve based 
on Newtonian dynamics would be mistaken as failing to explain the observed flat rotation curve when 
Keplerian dynamics were inappropriately employed, because it instead predicts a falling rotation curve. 



5 CONCLUSIONS 

In the present paper, we show that with appropriate mathematical treatments the apparent difficulties as- 
sociated with singularities in computing elliptic integrals can be eliminated when modeling Newtonian 
dynamics of thin-disk galactic rotation. Using the well-established boundary element techniques, the 
nondimensionalized governing equations for disks of finite sizes can be discretized, transformed into 
a linear algebra matrix equation, and solved by straightforward Gauss elimination in one step without 
further iterations. Although the mathematical derivations in Appendix A for removing the singularities 
seem somewhat sophisticated, the actual implementations of the numerical code are not as lengthy. With 
our code on a typical personal computer with a single Pentium 4 processor, each solution in Section 3 
takes no more than a minute or so to compute. Thus, a numerical code implemented according to our al- 
gorithm can be conveniently used to accurately determine the surface mass density distribution in a disk 
galaxy from a measured rotation curve (or vice versa), which is important for in-depth understanding 
of the Newtonian dynamics and its capability of explaining the "galaxy rotation problem" via rotation 
curve analysis. Moreover, the dimensionless galactic rotation number A in our mathematical system can 
provide important insights into the general dynamical behavior of disk galaxies. 

Through a systematic computational analysis, we have illustrated that the value of the galactic rota- 
tion number remains within ±10% of A = 1.70 for a wide variety of rotation curves. For most Sb type 
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galaxies like the Milky Way, having rotation curves typically with a very steep rise in a small central 
core region and a large range of flat portion, we have showed that A as 1.60 with a surface mass density 
very close to that of Mestel's disk (except in an infinitestmal neighborhood of the galactic center where 
the Mestel disk becomes singular). But for galaxies with "non-ideal" rotation curves containing consid- 
erable irregularities, our numerical approach can easily be used without modification for computing the 
corresponding surface mass density distributions accurately for rotation curve analysis. 

Because the value of A = Vq R g j (M g G) remains almost invariant for various galaxies, we can 
draw a conclusion that the total mass in a disk galaxy M g must be proportional to V 2 R g . For galaxies 
with similar characteristic rotation velocity Vb, their total mass M g must be proportional to their disk 
size R g . Our model predicts that at the disk edge the surface mass density is expected to diminish 
precipitously whereas within the disk edge the surface mass density should vary rather smoothly without 
sharp changes except near the galactic center. Thus, a disk galaxy with a finite amount of mass must 
also have a finite size, based on the Newtonian dynamics modeling. 

For a disk galaxy with a typical flat rotation curve, our modeling result show that the surface mass 
density monotonically decreases from the galactic center toward periphery, according to Newtonian 
dynamics. In a large portion of the galaxy, the surface mass density follows an approximately expo- 
nential law of decay with respect to the galactic radial coordinate. Yet the radial scale length for the 
exponential portion of surface mass density seems to be generally larger than that of the measured ex- 
ponential brightness distribution, suggesting an increasing mass-to-light ratio with the radial distance in 
a disk galaxy. This is consistent with typical edge-on views of disk galaxies often revealing a dark edge 
against a bright background bulge. 
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Appendix A: TREATMENTS OF SINGULAR ELEMENTS 

The com plete elliptic integrals of the fi rst kind and second kind can be numerically computed with the 
formulas (lAbramowitz & Stegunll 19721) 

4 4 

K(m) = ^a l m[-log(m 1 )^2b l m[ (A.l) 



and 



1=0 1=0 

4 4 



E(m) = 1 + 51 C;m i ~ lo s( TO i) dim i ' (A - 2) 

where 



i=i i=i 



(f — r\ 2 
• (A.3) 
r + r J 

Clearly, the terms associated with K(rrii) an d E(m,i) in (O become singular when f — > r, on the 
elements with r, t as one of their end points. 

The logarithmic singularity can be treated by converting the singular one-dimensional integrals into 
non-singular two-dimensional integrals by virtue of the identities: 



Jo /(£) log & = - Jo Jo f(£v)dridZ 

So f(0 i°g(i - £K = - Jo Jo /(! - 



(A.4) 
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where /(£) denotes a well-behaving (non-singular) function of £ on < £ < 1. 

But a more serious non-integrable singularity l/(f — n) exists due to the term E(rrii)/(r — r.;) 
in (O as f — > r,. The l/(r — ry) type of singu larity is treated by taking the Cauchy principle value to 
obtain meaningful evaluation (cf. | Kanwallll99^). as commonly done wi th the boundary element method 
dSladek & Sladek1ll998tl&ay1ll998tlSutradhar. Paulino & Grayll2008l) . In view of the fact that each n 
is considered to be shared by two adjacent elements covering the intervals [r,_i, r,] and [r^, r^+i], the 
Cauchy principle value of the integral over these two elements is given by 



lim 

r->() 



p(f)fdf r r *+ 1 p(f~jfdf 
f — r, 



r i+ i 

Ti-\-e 



(A.5) 



In terms of elemental £ 




is equivalent to 



l-e/(i-i-r»_i) 



[ft-i(l-Q+M]h-i(l-£) + r t £]d£ 
[Ml -0 + Pi+it][n(l - £) + rj+igdg 1 

Performing integration by parts on (1A.61 yields 

>i+i-M f /■ 1 d{[ /0i _ 1 (l-£)+p i £][r i _ 1 (l-£) + r i £]} 



(A.6) 



Pi n log 



d£ 



log(l-£)d£ 
log£d£ N 



where the two terms associated with log e cancel out each other, the terms with e log e become zero at the 
limit of e — > 0, and the first term becomes nonzero when the mesh nodes are not uniformly distributed 
(namely, the adjacent elements are not of the same segment size). In other words, inclusion of this first 
term enables the usage of nonuniformly distributed nodes for m ore effective computations, which is one 
of the algorithm improvements over that in our previous works dGallo & F eng 20091 1201 Oh . 
At the galaxy center r, = (i.e., i = 1), 



r i+ i 



p(f)df . 



(A.7) 



Thus, the l/(f — r.;) type of singularity disappears naturally. However, nume rical difficulty can still arise 
if p itself becomes singular as r — > 0, e.g., p cx 1/r as for the Mestel disk dMeste]||l963l) . The singular 
mass density at r = corresponds to a mathematical cusp, which usually indicates the need of finer 
resolution in the physical space. To avoid the cusp in mass density at the galactic center, we can impose 
a requirement of continuity of the derivative of p at the galaxy center r = 0. This be easily implemented 
at the first node i = 1 to demand dp/dr = at r = 0. In discretized form for n = Owe simply have 



p(n) = p(r 2 ) . 



(A.8) 



When ri = 1 (i.e., i = N), we are at the end node of the problem domain. Here we use a numer- 
ically relaxing boundary condition by considering an additional element beyond the domain boundary 
covering the interval [r», ri+i], because it is needed to obtain a meaningful Cauchy principle value. In 
doing so we can also assume r.; + i — r, = r, — r;_i such that log[(ri_|_i — ri)/(r, — becomes 
zero, to simplify the numerical implementation. Moreover, it is reasonable to assume pi + \ = because 
it is located outside the disk edge where the extremely low intergalactic mass density is expected to 
have inconsequential gravitational effect. With sufficiently fine local discretization, this extra element 
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can be considered to cover a diminishing physical space such that its existence becomes numerically 
inconsequential. Thus, at r; = 1 (where i = N) we have 



Now that only logarithmic singularities are left, (IA.4t can be used to eliminate all singularities in com- 
puting the integrals in ©. 

Noteworthy here is that the (removable) singularities in the kernels of the integral equation ©, 
when properly treated, lead to a diag onally dominant J acobian matrix with bounded condition number 
in the Newton-Raphson formulation JPress et al.lll988l) . This fact makes the matrix equation robust for 
any straightforward matrix solver. 
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